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Abstract 

Model reduction using the proper orthogonal decomposition (POD) method is 
applied to the dynamics of ferroelastic patches to study the first order square to 
rectangular phase transformations. Governing equations for the system dynamics 
are constructed by using the Landau-Ginzburg theory and are solved numerically 
By using the POD method, a set of empirical orthogonal basis functions is first 
constructed, then the system is projected onto the subspace spanned by a small set 
of basis functions determined by the associated singular values. The performance 
of the low dimensional model is verified by simulating nonlinear thermo-mechanical 
waves and square to rectangular transformations in a ferroelastic patch. Compari- 
son between numerical results obtained from the original PDE model and the low 
dimensional one is carried out. 

Key words: Phase transformation, ferroelastic patch, model reduction, proper 
orthogonal decomposition, Galerkin projection. 

1 Introduction 

The first order martensite transformation is a common transformation in many ferroelas- 
tic, ferromagnetic, and ferroelectric materials that holds the key to their specific prop- 
erties. For many new engineering applications of such materials, a better understanding 
of the mechanism of the transformation becomes decisively important. As a result, the 
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martensite transformations observed in smart materials have attracted substantial atten- 
tion from engineers, mathematicians, physicists, and control theorists. These materials 
possess unique properties of being able to sense, actuate, and inherently respond to exter- 
nal stimuli ([9] and references therein). Much attention is devoted to control and design 
optimization of these materials where mathematical modelling tools become ubiquitous. 
However, even for the systems governed by linear Partial Differential Equations (PDE), 
these issues are far from trivial [HI OH]- At the same time, mathematical models describing 
the dynamic behavior of most smart materials are intrinsically nonlinear, coupling differ- 
ent physical fields by conservation law equations and constitutive relations (e.g., [3ll5lfl2]). 
To bridge the gap between the nonlinear PDEs such as those arising in the description 
of the smart material systems and the low dimensional models suitable for control and 
optimization in realistic engineering applications, it is natural to attempt to approximate 
the dynamic behavior of PDEs by low dimensional models. Such low dimensional models 
can be obtained from the original systems by model reduction and the POD provides 
an efficient tool for doing so as soon as a collection of system states is available. This 
idea has been successfully applied to many active control problems involving fluid flows 
(see [21 IB] and references therein). Recently, we applied this idea to the analysis of phase 
transformations in the ID case [13J. In its essence, the POD methodology is close to the 
principal component analysis, techniques based on the singular value decomposition, and 
Karhunen-Loeve decomposition [21 IB]- 

In this paper, we develop a low dimensional model for the dynamics of ferroelastic 
patches exhibiting square to rectangular phase transformations to study 2D nonlinear 
thermomechanical waves. First, the original system of PDEs is analyzed numerically with 
Chebyshev's collocation method. Then, the basis functions are extracted from the ob- 
tained numerical results using the POD. Finally, the low dimensional model is constructed 
using the Galerkin projection method. The dynamic behavior of ferroelastic patches is 
simulated by the empirical low dimensional model, and its performance is evaluated by 
comparing the numerical results with those from the original PDE model. 

2 PDE Model and Its Numerical Solution 

Square to rectangular transformations are regarded as a 2D analog of the cubic to tetrag- 
onal or tetragonal to orthorhombic transformations observed in the general 3D case (see, 
e.g., [3] and reference therein). At mesoscale, there is a high temperature phase of greater 
symmetry (austenite, represented schematically as square), and two lower temperature 
phases (martensite variants, represented schematically as rectangles). The transforma- 
tion from austenite to martensite is not unique due to the fact that there are multiple 
equivalent orientations of the lower symmetry lattice with respect to the parent lattice. The 
transformation can be induced by either mechanical or thermal loadings (or their combi- 
nation), but only mechanically induced transformation will be modelled here. 

The governing equations describing the macroscopic dynamics of the ferroelastic struc- 
ture exhibiting the above mentioned transformation can be formulated on the basis of 
conservation laws. Following [12], we consider the following coupled system of PDEs 
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based on conservation of momentum and energy: 



d 2 u x _ dan + dcr 12 + , d 2 u 2 _ da 12 + da 22 + , 
^ dt 2 dx dy X1 ^ dt 2 dx dy y ' 

de , fd 2 e d 2 e\ n de 2 
Cv m =k W + df) +a2ee2 ^ + g ' 

a n = -^{axe x + a 2 {6 - 9 )e 2 - a A e\ + a 6 e|) + y V^e 2 , W 
1 

0"l2 = = CT 2 1, 

o"22 = -^{axei - a 2 (6 - 9 )e 2 + a 4 e\ - a e e 5 2 + y Vje 2 ), 

where c v is the specific heat constant, 6q is the reference temperature for the trans- 
formation, V is the gradient operator, a\ a 2 , a 3 , a^, a 6 , and d 2 are the material-specific 
coefficients, and e±, e 2 , e 3 are dilatational, deviatoric, and shear components of the strains, 
respectively. f x , f y and g are mechanical and thermal loadings, respectively. 

In the above model, e 2 is chosen as the only order parameter to characterize differences 
in phases. Strain components are defined as follows: 



ei = (7711 + 7722) /v2, 

e 2 = (7711 - 7722) / ^ ( 2 ) 
e3 = O/12 + V21) /2. 

where the Cauchy-Lagrangian strain tensor rj is given by its components as follows(with 
the repeated- index convention used): 

j(y(xi(H ^ + ^ A (3) 

In the above formulation, stress components are kept as independent variables to be 
solved for. This approach is similar to the DAE approach first proposed in the context of 
shape memory alloys materials modelling in 



2.1 Chebyshev Collocation Methodology 

The system formulated in the previous section is solved here numerically by the method of 
lines. First, we discretize the system spatially, converting it into a system of differential- 
algebraic equations, and then we apply a backward differential formula to integrate the 
resulting system in time. 

For the spatial discretization, the Chebyshev pseudospectral approximation is em- 
ployed on a set of 2D Chebyshev points (xi,yj) G [—1, 1] x [—1, 1] pQ: 

Xi = cos(— ),yj = cos(— ), i,j = 0, l,...,iV, (4) 
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where N+l is the number of nodes. Using the discretization, all the unknown distributions 
in the patch can be approximated by the following linear combination: 

N N 

f(x, y) = Yl Mi( x )My)> ( 5 ) 

where f(x,y) could be either stresses, displacements, velocities, or temperature; f it j is 
the function value at (xi,yj); 4>i(x) and 4>j(y) are the i th and j th Lagrange interpolating 
polynomials in x and y directions, respectively. 

Because Eq.flSJ) is just a linear combination of the interpolating polynomials, the deriva- 
tives of the functions df(x,y)/dx and df(x,y)/dy can be easily obtained by calculating 
d(f>i(x)/dx and d<fij(y)/dy. Following the idea given in[TT], all the differential operators 
in Eq.flT]) can be written SIS db matrix form: 



d f( x ,y) n jp d f( x ,y) 

= U T t , 

ox 



= DyF, (6) 

Vi 



xi,yj dy 

and similarly for the higher order differentiation operators. The differentiation matrices 
D x , D y , D xx , and D yy can be calculated using the approximation given by Eq.(j3]). 

By approximating all the differential operators in the governing Eq.flT]), the system 
can be discretized and converted into a set of nonlinear ODEs and algebraic equations. 
The boundary conditions of the system can be easily enforced [TT] [1] . The system after 
discretization finally takes the following form: 

M^ + G(t,X,U) = 0, (7) 

where X is a 5(iV + 1) vector collecting all the unknowns we are seeking for, the matrix 
M = diag(ai, d2, 05(jv+i)) is a 5(iV + 1) x 5(N + 1) matrix having entries "one" for all 
the differential equations associated with ODEs and "zero" for all algebraic equations. 

The above DAEs system is a stiff system and must be solved by an implicit algorithm. 
Here, the second order backward differentiation formula method is employed for this 
purpose. By discretizing the time derivative using the second order approximation, the 
DAE system can be converted into the following algebraic system at each time level: 

M 0X n - 2X- 1 + ^X n - 2 ) + AtN (t n , X n , V(t n )) = 0, (8) 

where n denotes the current computational time layer. For each computational time layer, 
iterations are carried out with the bi-conjugated gradient iteration method for computing 
X n having X n_1 and X n ~ 2 . Starting from the initial value, the vector of unknowns X is 
sought for all specified time instances employing this algorithm. 

We note that to deal with strong (cubic and quintic) nonlinearities in the order pa- 
rameter, a smoothing procedure similar to that proposed in [6] has been employed during 
the iteration process. In particular, we have used the following expansions: 

3 5 

y'^E^ ^gE^-l: (9) 

1=0 i=0 

where y stands for €2 here, and n time layer. 
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3 Construction of Low Dimensional Models 



3.1 Proper Orthogonal Decomposition 

One of the most effective ways to construct a low dimensional model for approximation 
to a dynamical system is to find a set of optimal basis functions by which the dynamical 
system can be approximated best at any time. In what follows, we apply the POD 
method to construct an empirical basis for the system dynamics describing 2D phase 
transformations. 

For the dynamical system given by Eq.(TjQ), the system state at any time can be fully 
characterized by the displacements, velocity, and temperature distributions, while dis- 
placements and velocities can be represented by using the same basis functions. Here, we 
only need to construct spatial basis functions for displacement and temperature distribu- 
tions. Let us assume that U(x,t), which describes the system dynamics, belongs to an 
inner product (Hilbert) space 7i. The POD is concerned with the possibility to find a set 
of orthonormal spatial basis functions (f) = {<f)j(x),j = 1, . . . , P} in 7i which are optimal 
in the sense that the P dimensional approximation 

p 

U P {x,t) = Y,<t)<p i {x) (10) 

i=l 

gives the best (in the least square sense) approximation to the function 14 (x, t) among all 
those P dimensional approximations in H, [21 E]- Here are the general Fourier coefficients 
associated with {<fii}. They are functions of time. Mathematically, the idea of POD is to 
choose the basis functions {<fik} to maximize the mean projection of the function U(x,t) 
on {<fr k } 

- ("> 

<pcL 2 (n) \\(j)\\ 2 

where E(-) denotes the mean value functional, | • | is the modulus, ( • ) is the inner 
product, and || • || is the L 2 norm [21 E]- 

For convenience of the projection operation, that will be employed later, the orthonor- 
mality requirement is enforced on the basis functions we are looking for: 

<*.*>={ i li = A (12) 

The general Fourier coefficients in this case can be calculated simply as: 

a i (t) = (U(x,t),4 >i ). (13) 

Taking into account the above requirement on the basis functions, the maximization 
problem can be reformulated in terms of calculus of variation, with a functional for the 
constrained variational problem [21 E] : 

J[0] = E(|(W,0)| 2 )-A(||0|| 2 -1). (14) 
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A necessary condition for extrema is that the first order variation of the functional with 
respect to all variations <f> + 8i[) should be zero: 



^J[0 + <fy]|«=o = O, 
which finally leads to the following eigenvalue problem (2J [7j : 



(15) 



/ E(U{x)U(x'))<p{x')dx' = \<j>{x) 



(16) 



Jn 



whose kernel /C = E(U(x)U(x')) is the auto-covariance function of the two points x and 
x'. 

For the current problem, the dynamics of the system can be sampled at any time 
instances we need, by using the numerical algorithm mentioned above for Eq.flTJ. Let us 
denote by U l the system state at the i th time instance, called the i th snapshot. In discrete 
form, each snapshot can be written as a column vector with M entries, where M is the 
number of nodes for spatial discretization. 

To construct the orthonormal basis, all the snapshots are collected in one matrix 
U = {If 1 ,i = 1, . . . , N}. If one put each snapshot as one column in the collection matrix, 
the collection U will be a matrix of M x N. The orthonormal basis vectors for the given 
collection U can be calculated by the singular value decomposition as follows: 



where L is M x M orthonormal matrix, R is a N x N orthonormal matrix and the 
superscripts T indicates matrix transpose. S is a M x N matrix with all elements zero 
except along the diagonal, those non-zero elements are arranged in a decreasing way along 
the diagonal. They are the singular values of U, emphasizing the relative importance of 
its associated eigenvectors in L and R. 

Let SR T = Q T in the singular value decomposition, then U = LQ T . Let 4>k be the k th 
column of L and be the k th row of Q, then the matrix U's singular value decomposition 
can be rewritten as 



where m = min(N, M) is the rank of the collection matrix U. This approximation is what 
we are looking for in Eq.f llOl) . The lower dimensional approximation of the matrix U can 
be easily obtained by keeping the first few largest singular values and their associated 
eigenvectors in L and Q. The number of eigenvectors should be determined by compro- 
mising between the dimension number of the resultant system and the approximation 
accuracy [3, 12] • 

3.2 Galerkin Projection 

It has already been mentioned that Eq.flT]) can be converted into a system of DAEs like 
Eq.(j7|) by discretization. It can also be converted into a system of ODEs if we substitute 



U = LSR T 



(17) 



ni 




(18) 



k=l 
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all the constitutive relations into Eq.flT]). The latter approach is applied here for the 
construction of the low dimensional model. For convenience, Eq.([T]) is converted into a 
dynamical system representation as follows: 

j t X = Q(t,X,H), (19) 

where Q is the collection of functions given by the spatial discretization of Eq.Q, while 
H is the collection of inputs after conversion. Vector x collects all unknowns includ- 
ing displacements in the x and y directions, velocities in the x and y directions, and 
temperatures in the collocation nodes. 

Having obtained the optimal basis functions for the dynamical system, its lower di- 
mensional approximation can be obtained by projecting the full system orthogonally onto 
the subspace spanned by the chosen basis functions. It is known that the model reduction 
can be achieved due to the fact that much smaller number of basis functions are needed 
to approximate the full system, when the chosen basis functions are optimal. 

Let us denote the projection operator as P r and Pj its transpose. Then, a low di- 
mensional system on the chosen subspace S can be constructed by using the following 
rule: 

For any point Z G S, compute the vector field Q(t, H , Pj ' Z) and take the projection 
P r Q(t, H, Pj Z) on to the subspace S. It should be equal to Z. The result can be 
formulated as: 

j t Z = P r Q{t,PjZ,H). (20) 

The approximation will introduce an error defined as follows: 

r = X - PjZ. (21) 

To achieve the best approximation with the given basis, the error function of the 
approximation is required to be orthogonal to all the basis functions at any time: 

(r,<p k )= / r(x)(f)(x)dx = 0. (22) 
Jn 

It is easy to check that the low dimensional system given by Eq.( l20l) satisfies this condition, 
if the orthogonality of the basis functions is used. 



4 Numerical Results 

To demonstrate the performance of the proposed methodology, numerical results from 
the low dimensional model for the dynamics of a ferroelastic patch are compared with 
its counterparts from the PDE model. The dynamical behavior of a Au23Cu3oZn 47 patch 
with size of [0, 1] x [0, l]cm 2 is analyzed numerically. All the physical parameters for this 
specific material are available in the literature (e.g., [U [12]). It has been shown that 
the first order square to rectangular transformation can be mechanically induced in this 
material [U EJ H2]- It has also been shown that the nonlinearity is very strong, and the 
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nonlinear coupling between the mechanical and thermal fields complicates the problem 
further. 

The strategy here is as follows: 
We first perform the numerical simulation using Eq. (JTJ) with a representative mechanical 
load. As a result, the collection matrix can be constructed and the empirical eigenfunc- 
tions can be extracted. Then, we simulate the dynamics using Eq. f[2TJ]) with loads that 
differ from the previous loading. By comparing the numerical results from Eq. f[2"U]) with 
those from Eq.flT]), we assess the performance of the low dimensional model. 

The initial conditions for all the simulations reported below are set the same. They 
are 9(x, 0) = 250° K and u\ = u% = v\ = t> 2 = 0. Boundary conditions for Eq.(TjQ) are 
taken as pinned-end mechanically and insulated thermally. The distributed mechanical 
loading for collecting snapshots is the following (for one period, gj (ms 2 cm 2 )) 

0<t<2 
2 < t < 4 
4 < t < 6 
6 < t < 8 
8 < t < 10 
10 < t < 12. 

There are 13 nodes, each used for «i, u 2 , t>i, t>2, and 9 discretization. Two periods of 
loadings are performed (t G [0,24]ms) and totally 201 evenly distributed snapshots are 
sampled. The collection matrix for each field variable then can be constructed by collecting 
their values on each sampling instance. To demonstrate the phase transformation, the 
order parameter and temperature on the central horizontal line are plotted versus time 
in FigHJ 

With the above constructed collection matrix at hand, the POD for each of the field 
variables can be done by using Eq.t flTj) and Eq.(fl8l). by which the eigenfunctions for 
Mi, w 2 , f i, t>2 an d 9 can be estimated [13]. The first three eigenfunctions for displacement 
Ui are plotted in the left column of FigfSl together with those three for temperature in the 
right column. The eigenfunctions for v i can be taken as the same as those for u\ because 
their differences are in time domain, and similar for t>2. 

For the construction of the low dimensional model, the first 12 basis functions for 
Mi,M2, and 9 are kept, so the system given by Eq. (1201) has a dimension of 60, 12 for each 
of the variables ui, 112, fi, t>2 and 9 each. A standard ODE integrator (ode45 in Matlab) 
is applied to simulate the state evolution. 

The mechanical loadings now are f x = f y = 6000 sin(7rt/6), and the simulation is 
performed within the same time interval. The simulated distributions of e<i and 9 from the 
low dimensional model are presented in the top row of Figj3j The simulated distributions 
clearly indicate the mechanically induced phase transformation, and the temperature 
oscillation driven by the mechanical loads due to the thermomechanical coupling. To 
show the effect of transformation, the distribution of €2 at t = 3 and t = 9 from the low 
dimensional model is presented in the middle row of the figure. In this case, the entire 
patch is divided into two sub-areas, associated with the two martensite variants. The 



7000t/2 
7000(4 - 1)/2 


' 7000(6 -t)/2 
7000(t - 10)/2 
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combination of the two martensite variants also depends on the loading direction [12] . To 
demonstrate the thermo-mechanical coupling, the temperature distributions 9 at t = 3 
and t = 9 also plotted in the figure (the bottom row). It can be observed that after the 
transformation, the temperature distribution of the patch is also changing. 

To validate the low dimensional model, a more detailed comparison between the two 
sets of numerical results (from the PDE model and POD model ) is carried out here. 
To do so, three points in the shape memory alloy patch are chosen, and the strain and 
temperature values from the two models are compared at these points. The first point is 
chosen close to the left boundary as (0.067,0.5), the strain (left) and temperature (right) 
evolutions at the point during the loading process are presented in the top row of FigJU 
The second point is chosen close to the center of the patch as (0.5,0.63), and the strain 
and temperature evolutions are presented similarly in the middle row of the figure, while 
those for the third point (0.75, 0.25) are presented in the bottom row. At all the three 
points, the strain evolutions simulated using the low dimensional model agree very well 
with those from the PDE model, the discrepancies can be distinguished only by a close 
look at the plots. The phase transition is indicated clearly by the strain evolution. 

The temperature evolutions from the low dimensional model at the three points also 
agree qualitatively with their counterparts from the PDE model, the oscillatory behaviours 
of the temperature driven by the mechanical loadings due to the thermo-mechanical cou- 
pling are also successfully captured in the low dimensional model. The agreement is not 
perfect in temperature evolutions, there are certain discrepancies between the two mod- 
els. This can be explained as follows. In the dynamics of the patch, the thermal and 
mechanical fields are coupled in the PDE model, but for the low dimensional model it 
is not practical to extract the eigenmodes for the thermal and mechanical fields simul- 
taneously by using the collection matrix including both temperature and displacement 
distributions, because their magnitude is very different ( 250 versus 0.1), and character- 
istics of the variable with smaller magnitude will be ignored if one do so. In the current 
paper, the eigenmodes for the displacement distributions are obtained by the POD anal- 
ysis using the collection matrix of displacements only, and similarly for the eigenmodes of 
temperature distributions. With this simplification, some characteristics associated with 
the coupling between the thermal and mechanical fields are lost. This loss will induce 
minor errors in the low dimensional approximation. When the dynamics of the patch are 
simulated using the low dimensional model with a mechanical loading, the characteristics 
of the thermal behaviour driven by the mechanical loading can be captured, but the minor 
error due to the partly loss of coupling characteristics will be accumulated and make the 
temperature evolutions slightly different from those of the PDE model, just as indicated 
by the comparison in FigJU 

From the above numerical experiments and comparison, it is clearly shown that the 
low dimensional model is able to capture the characteristics of the dynamics involving 
phase transformation, as well as the thermo-mechanical coupling. 
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5 Conclusion and Discussion 



In this paper, the dynamics of a ferroelastic patch exhibiting square to rectangular marten- 
site transformations has been analyzed. The system of PDEs based on the Landau theory 
of phase transition has been analyzed numerically by the Chebyshev collocation method. 
The proper orthogonal decomposition has been employed to extract the eigenmodes of 
the dynamics of the patch. A low dimension system has been constructed by projecting 
the PDEs system onto the subspace spanned by the first few eigenmodes. 

Numerical experiments demonstrated that the complicated dynamics of ferroelastic 
materials can be modelled successfully with the derived low dimensional system. 

References 

[I] Q. Alfio, S. Riccardo, and S. Fausto, Numerical mathematics, Springer- Verlag, Berlin, 
2000. 

[2] P. Homes, J. L. Lumlay, G. Berkooz, J. C. Mattingly, and R. W. Winttenburg, Low di- 
mensional structure in coherent structures in turbulence, Physical Reports 287 (1997) 
337-384. 

[3] T. Lookman, S. R. Shenoy, K. O. Rasmusseh, A. Saxena, and A. R. Bishop, Ferroe- 
lastic dynamics and strain compatibility, Physical Review, B67 (2003) 024114. 

[4] R. Melnik, A. Roberts, K. Thomas, Phase transitions in shape memory alloys with hy- 
perbolic heat conduction and differential algebraic models, Computational Mechanics 
29 (1) (2002) 16-26. 

[5] R. Melnik, A. Roberts, K. Thomas, Coupled thermomechanical dynamics of phase 
transitions in shape memory alloys and related hysteresis phenomena, Mechanics Re- 
search Communications 28 6 (2001) 637-651. 

[6] M. Niezgodka, J. Sprekels, Convergent numerical approximations of the thermome- 
chanical phase transitions in shape memory alloys, Numerische Mathematik, 58 (1991) 
759-778. 

[7] M. Rathinam, L. R. Petzold, A new look at proper orthogonal decomposition, 54AM, 
J. Numer. Anal. Vol 41 5 (2003) 1893-1925. 

[8] C. W. Rowley, T. Colonius, R. M. Murray, Model Reduction for Compressible Flows 
Using POD and Galerkin Projection, Physic D 189 (2004) 115-129. 

[9] D. K. Sood, R. A. Lawes, V. V. Varadan (Eds), Smart Structures and Devices, Pro- 
ceedings of SPIE, December 2000, Melbourne, Australia. 

[10] R. Teman, Infinite Dimensional Dynamical Systems in Mechanics and Physics. 
Springer- Verlag, New York, 1998. 

[II] L. N. Trefethen, Spectral method in Matlab, SIAM Philadelphia, 2000. 



10 



[12] L. X. Wang, and R.V.N. Melnik, Thermomechanical waves in SMA patches under 
small mechanical loadings, in Lecture Notes in Computer Science 3039, (2004), 645- 
652. 

[13] L. X. Wang, and R.V.N. Melnik, Simulation of thermomechanical wave with an empir- 
ical low dimensional model, Lecture Notes in Computer Science 3514, (2005), 884-891. 



11 



0.2 




Figure 1: Simulated evolution of strain and temperature on the central horizontal line 
in the patch involving phase transformation, using PDE model, (top) Strain evolution, 
(bottom) Temperature evolution. 



12 



Mode 1 



Mode 1 




POD model 



POD model 




Figure 3: Numerical simulation of the dynamics of a shape memory alloy patch using a 
low dimensional model. From top to bottom, left to right, (a) Strain evolution on the 
central horizontal line, (b) Temperature evolution on the central horizontal line, (c) Strain 
distribution at t=3, (d) Strain distribution at t=9, (e) Temperature distribution at t=3, 
(f) Temperature distribution at t=9. 
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Time Time 

Figure 4: Comparision of numerical results from the low dimensional model with those 
from the PDE model. From top to bottom, left to right, (a) Strain evolution at (0.067,0.5), 
(b) Temperature evolution at (0.067,0.5), (c) Strain evolution at (0.5,0.63), (d) Tempera- 
ture evolution at (0.5,0.63), (e) Strain evolution at (0.75,0.25), (f) Temperature evolution 
at (0.75,0.25). 
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